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Tidal mass loss in star clusters and treatment of escapers in 
Fokker-Planck models 
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ABSTRACT 

This paper presents a new scheme to treat escaping stars in the orbit-averaged Fokker-Planck 
models of globular star clusters in a galactic tidal field. The existence of a large number of 
potential escapers, which have energies above the escape energy but are still within the tidal 
radius, is taken into account in the models. The models allow potential escapers to experience 
gravitational scatterings before they leave clusters and thus some of them may lose enough 
energy to be bound again. It is shown that the mass evolution of the Fokker-Planck models 
are in good agreement with that of TV-body models including the full tidal-force field. The 
mass-loss time does not simply scale with the relaxation time due to the existence of potential 
escapers; it increases with the number of stars more slowly than the relaxation time, though it 
tends to be proportional to the relaxation time in the limit of a weak tidal field. The Fokker- 
Planck models include two parameters, the coefficient 7 in the Coulomb logarithm ln(7iV) 
and the coefficient v c controlling the efficiency of the mass loss. The values of these param- 
eters are determined by comparing the Fokker-Planck models with the iV-body models. It is 
found that the parameter set (7, v c ) = (0.11, 7) works well for both single-mass and multi- 
mass clusters, but that the parameter set (7, v c ) = (0.02, 40) is another possible choice for 
multi-mass clusters. 
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1 INTRODUCTION 

The numerical integration schem e of the orbi t-averaged Fokker- 
Planck (FP) equation developed bv lCohnl i ll 9791) has been one of the 
most useful tools for simulating the dynamical evolution of globu- 
lar star clusters. In addition to two-body relaxation, many physical 
processes have been incorporated into FP models to achieve real- 
istic modelling of the globular cluster evolution; these processes 
include tidal cutoff, binary heati ng, disc and bulge shocks, m ass 
loss via stellar evolution, etc. (see lShin. Kim & Ta kahashi 2008 for 
a recent example of detailed FP modelling). 

In this paper we consider the dynamical evolution of globu- 
lar clusters in a steady galactic tidal field. Our main purpose is to 
investigate what boundary condition can give a better description 
of escape of stars from clusters in the tidal field. Thi s stud y has 
been motivated b y the studies bv iFukushige & Heggiel J20001) and 
iBaumgardtl J200lb. 

iFukushige & Heggiel feOOOh found that a large fraction of stars 
with energies above the escape energy (i.e. potential escapers) take 
much longer escape time than the dynamical time. Until their study 
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it had been generally thought that the escape time-scale is of the or- 
der of the dynamical time and that the mass-loss times of the clus- 
ters essentially scale with the relaxation time, which is much longer 
than t he dynamical time. The findings of Fukushige & Heggie 
(2000) indicate that this simple scaling may be spoiled by poten- 
tial escapers with long escape times. 



In fact[Baumgardt| feOOll) performed iV-body simulations and 
showed that the mass-loss times (lifetimes) of clusters do not scale 
with the relaxation time i r h but scale with t^h • He concluded 
that the reason is that some of potential escapers are scattered 
back to lower energies before they leave the cluster. More recently 
Tanik awa & Fukushigd J200II) showed that the dependence on the 
relaxation time changes with the strength of the tidal filed. These 
two studies have revealed that the behavior of potential escapers 
greatly influences the rate of mass loss from clusters in the tidal 
field. 

The effects of long escape times and re-scattering of potential 
escapers have never been considered in previous FP models in the 
literature, but it was assumed that escapers leave a cluster on the 
dynamical time-scale, as is described in detail in Section [2] Since 
the effect of the galactic tidal field is essentially important to the 
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cluster evolution, it is necessary to find a way to include the effect 
into FP models as precisely as possible. 

We should mention that Takahashi & Portegies Zwart (1998, 
2000) compared FP and TV-body models of star clusters in the 
tidal field and found good agreement between these two theoret- 
ical models over a wide range of initial conditions. They showed 
that the use of anisotropic FP models w ith the apocentre escape 
criterion dTakahashi. L ee & Inagaki 1997) and the dynamical-time 
removal of escapers dLee & Ostrikenll987l) is necessary to obtain 
such good agreement. However, note that in their TV-body models 
the tidal force field is not included but the tidal cutoff is applied. 
Takahashi & Portegies Zwart (2000) confirmed that the difference 
between tidal cutoff and self-consistent tidal field TV-body models 
is small for a particular set of initial conditions, but did not do sys- 
tematic investigations on this problem. 

In this study we have devised a new scheme to treat escapers 
in FP models. The scheme defines a region of potential escapers 
in phase space and allows them to be scattered again. Comparing 
the results of FP models calculated with the new scheme with the 
results of TV-body models, we examine the accuracy of the FP mod- 
els. 



2 FOKKER-PLANCK MODELS OF STAR CLUSTERS IN 
A STEADY TIDAL FIELD 

2.1 Basic assumptions 

The orbit-averaged FP equation is de rived under the assumption 
of spherical symmetry of star clusters dCohnll 19791) . Therefore the 
tidal field, which is not spherically symmetric, cannot be directly 
incorporated into orbit-averaged FP models. In FP models the ef- 
fect of the tidal field is taken into account by imposing a tidal cutoff 
radius r t on the cluster, which is treated as an isolated system in 
other respects. Under these assumptions the distribution function / 
of stars at time t depends only on the energy of a star per unit mass, 
E, and the angular momentum per unit mass, J. 



2.2 Classical treatments of escapers 

First we summarise classical treatments of escapers used in FP 
models of previous studies. 



2.2.1 Escape criteria in phase space 

In previous studies, two kinds of criteria were adopted to define an 
escape region in (E, J)-space: 

(i) Energy criterion 

GM 



E > E t 



n 



(ii) Apocentre criterion 

r a (-E, J) > r t , 



(1) 



(2) 



where M is the cluster mass and r a (E, J) is the apocentre radius 
of a star having energy E and angular momentum J. It is assumed 
that a star is destined to escap e once it enters into the escape region. 

The apocentre criterion dTakahashi et alJI 19971) is considered 
to be more realistic, at least as long as the tidal field is mod- 
elled as a radial cut-off, and in fact gives better agreement be- 
tween FP and TV-body models (Takahashi & Port egies ZwarJ 19981 
2000). For isotropic FP models, where the distribution function 



does not depend on J , only the energy criterion can be applied (e.g. 
iLee & Ostrikeill 19871) . 



2.2.2 Removal of escapers 

In previous studies stars in the escape region are assumed to leave 
the cluster inevitably, as mentioned above. It is also assumed that 
the time required for this travel is of the order o f the dynamical 
time a t the tidal radius. Considering this travel time. lLee & Ostrikeil 
dl987l) applied the following equation to the distribution function / 
in the escape region: 



dt 



(I)' 



1/2 



/ttid 



(3) 



where v c is a di mensionless constant determ ining the efficiency of 
escape (see also lLee, Fahlma n & Richer 1991). The time-scale t t id 
is an orbital time-scale at the tidal radius defined by 



ttid 



2tt 



(4) 



where p t is the mean mass density within the tidal radius. 

Since the dynamical time is generally much smaller than the 
relaxation time in globular clusters, we may assume that escapers 
leave the cluster immediately after they enter into the escape region, 
when we are interested only in the evolution on the relaxation time- 
scale. This assumption leads to the boundary condition 

/ = (5) 
on the tidal boundary (e.g. IChernoff & Weinberglll99d) . 

2.3 A new treatment of escapers 

The boundary condition of equation ([3} takes account of the fact 
that stars satisfying the escape criterion, i.e. potential escapers, 
need time to actually leave the cluster. However the effect of re- 
scattering of potential escapers is not considered there. Here we 
propose a new scheme in which the re-scattering effect is taken 
into account. 

First we summarise basic assumptions and equations. Suppose 
that the cluster is on a circular orbit, with radius Rq and angular 
velocity cj, round the centre of a spherical galaxy. We consider the 
motion of a star in the rotating coordinate system moving with the 
cluster; the origin is at the cluster centre, the a>axis points to the 
galactic centre, and the y-axis is in the cluster orbital plane. If the 
cluster and the galaxy are treated as point masses M and Mg (S> 
M) and the size of the cluster is much smaller than Rq, there exists 
a conserved quantity known as the Jacobi integral given by 



E 3 



1 2/ 2 2n 

-uj (6x — z ) 



(6) 



/ _ GM 
2 ~ 2 

(cf. lSpitzerl[l987L Chapt. 5). Here v is the velocity of the star mea- 
sured in the rotating frame, r is the distance from the star to the 
cluster centre, and the angular velocity uj is given by 



GM G 



(7) 



The third term on the right-side in equation {6} is a combination of 
the centrifugal and tidal potentials. 

The effective potential is defined as 



x GM 1 a a 

es{x,y, z) = -uj {3x 

r 2 



(8) 
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A contour plot of <f> e g is shown, e.g., in Fig. 5.1 o f lSpitzeil fl987h . 
The effective potential has the saddle points at (±x e , 0, 0), where 



= (JL) 1/3 Rg 

i / 



V3M G , 



and 



4> ca (±x c , o,o) = — 



3 GM 

2 x e 



(9) 



(10) 



The equipotential surface passing through these saddle points in- 
tersects with the y-axis at y — ±y e , where 



Ue — jj^c 



(ID 



The necessary condition for escape of a star from the cluster is 
given by 

3 GM 
= ~2~' 



Ej > E 3 



(12) 



Note that equations d 1 0b . i ll It . and d!2b are valid for any spherical 
galactic potential. Fukushige & Heggie] J2000h found that the time- 
scale for escape of stars with Ej > Ej lC rit varies as 

-2 



t c OC - £j,crit)" 



(13) 



With this relation in mind we have devised a new scheme to 
follow the evolution of potential escapers. In this scheme the evolu- 
tion of the distribution function / for potential escapers is described 
by 

d.f fdf\ f 



df = (d£\ 
dt V dt ) , 



(14) 



' coll t a (E) 

where the first term on the right-side is the FP collision term and 
the second term represents mass loss due to escape. Here the escape 
time-scale t c is given by 

t e (E) t t id V -Ecrit / ' 

where v Q is a dimensionless numerical constant. It should be noted 
that energy E, not the Jacobi integral Ej, is used in equations dl4b 
and d!5t , Energy E does not include the centrifugal and tidal po- 
tentials. Despite this difference, we use the same critical value of 
energy 

„ 3 GM 

E CT it — — , 

2 r t 

where the tidal radius rt is identified with :r c . One might think that 
using equation d!5b with equation d!6b is too crude an approxima- 
tion, but it brings good agreement between FP and TV-body models 
as is shown in Section[5] 

The most important difference between equations {3j an d {33} 
is that the latter includes the collision term. Thus equation d 14b al- 
lows potential escapers to be scattered back to lower energies. The 
effect of mass loss is included in both equations in a similar way, 
though the functional forms of the escape time-scale t c are differ- 
ent. 

In this new treatment of the tidal field, the escape criteria de- 
scribed in Section [2.2.1l are modified as follows: 



(16) 



(i) Energy criterion 

E > E cr it = — — , 

2 n 

(ii) Apocentre criterion 
r*(E,J) > H rt . 



(17) 



(18) 



Note that eff (O, ±2r t /3, 0) = 0(0, ±2r t /3, 0) = -3GM/2r t . 
Equation J 1 4b is applied only in the region where an adopted crite- 
rion is satisfied. 



2.4 The Fokker-Planck code 

The FP code used in the present study is e ssenti ally the same as 
that used by Takahashi & Portegies Zwart (2000), but adopts the 
new scheme for treating escapers described above. The code cal- 
culates the evolution of the di s tributi on function f(E, J, t). Unlike 
lTakahashi~& Portegies Zwart] d2000t) . stellar evolution is not con- 
sidered in the models presented in this paper. Instead the effect 
of heating by three-body b inaries is considered in the manner de- 
scribed in Takahash] dl997h . 

For all the models presented in the present paper, 201 energy 
mesh points, 51 angular-momentum mesh points, and 101 radial 
mesh poin t s are u sed. The meshes are constructed as described in 
iTakahashil dl995h . When calculating the evolution of multi-mass 
clusters, 10 discrete mass-components are used to represent a con- 
tinuous mass function. 

Our FP models have two free parameters: one is v c in equa- 
tion d!5b and the other is 7 in the Coulomb logarithm ln(7TV) 
appearing in the FP collision term. How the val ue of v e is deter- 
mined is described in Section[3] We set 7 = 0.11 ( Giersz & Heggie 
Il994al) in most of our runs and 7 = 0.02 {g iersz & Heggie 1996) 
in a part of runs for multi-mass clusters. 



3 RESULTS 

3.1 Comparison with TV-body models: single-mass clusters 

First we compare FP models with the full tidal field models of 
iBaumgardtl j200ll) and additional TV-body runs performed for this 
comparison. All the model clusters are composed of equal-mass 
stars and move on circular orbits round a point-mass galaxy. The 
initial distribution of stars is given by King models jKingj[T 966). 

Results are presented in iV-body units, where the initial total 
mass and energy of a cluster are equal to 1 and —0.25, respectively, 
and the gravitational constant G = 1. The same units are used 
throughout this paper. 

Here we will refer to FP models with the boundary condition 
of equation J 14b as "FPf" models, which aim to model clusters in a 
self-consistent full tidal field. FP models with equation ((3) will be 
called "FPd" models, where stars beyond the tidal cutoff radius are 
removed on the dynamical time-scale. 

Fig.Q]compares FPf and TV-body models concerning the evo- 
lution of the total mass of bound stars. The initial models are 
Wo = 3 King models with the number of stars N = 1024, 
4096, 16384 and 65536. The new treatment of escapers described 
by equation U4t with the apocentre criterion of equation d 1 8b is 
employed in the FPf models. The agreement between the FPf and 
iV-body models is good in all the cases. In fact the value of the 
parameter u c in equation f 14b has been determined so that good 
agreement is obtained b y performing test runs with differe nt val- 
ues of v c as was done bv lTakahashi & Portegies Zwartl feOOd) . We 
have finally chosen the value of u c = 7. All the FPf models shown 
in Fig.[T]are calculated with this value. 

Fig. H] shows the evolution of the ratio of the mass of poten- 
tial escapers M pe to the total cluster mass M for the runs shown 
in Fig. Q] The agreement between the FPf and TV-body models is 
fairly good also in this comparison. Note that here M pc for the FPf 
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Figure 1. Evolution of the cluster mass. The solid lines represent FPf mod- 
els, and the dashed lines represent TV-body models. The initial models are 
W = 3 King models with the number of stars TV = 1024, 4096, 16384 
and 65536. 
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Figure 2. Evolution of the ratio of the mass of potential escapers M pc to 
the total cluster mass M. The ratio is plotted as a function of the cluster 
mass at each instance. 
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Figure 3. Half-mass time T^aif as a function of the initial half-mass relax- 
ation time t r h,i for the initial conditions of Wo = 3 King models. Two 
types of FP models, FPf and FPd models (see text), are shown by the circles 
and crosses, respectively, and TV-body models are shown by the triangles. 
The dotted lines represent scalings proportional to ar, d tfht ( me Y ale 
arbitrarily shifted in a vertical direction). 



Table 1. Half-mass times T^aif given by TV-body, FPf, and FPd models for 
the initial conditions of King models with Wo = 3. 



TV 


'rh. 




Thalf 


Thalf 


^half 










(TV-body) 


(FPf) 


(FPd) 


128 


5.13 


X 


10° 


8.94 x 10 1 


8.37 


X 


10 1 


6.76 x 10 1 


256 


8.13 


X 


10" 


1.27 x 10 2 


1.21 


X 


If) 2 


9.93 x 10 1 


512 


1.35 


X 


10 1 


1.83 x 10 2 


1.71 


X 


10 2 


1.50 x 10 2 


1024 


2.30 


X 


10 1 


2.59 x 10 2 


2.51 


X 


10 2 


2.41 x 10 2 


2048 


4.01 


X 


10 1 


3.73 x 10 2 


3.73 


X 


10 2 


4.00 x 10 2 


4096 


7.11 


X 


10 1 


5.58 x 10 2 


5.56 


X 


10 2 


6.86 x 10 2 


8192 


1.28 


X 


10- 


8.41 x 10 2 


8.33 


X 


10 2 


1.20 x 10 3 


16384 


2.32 


X 


10 2 


1.18 x 10 3 


1.26 


X 


10 3 


2.16 x 10 a 


32768 


4.24 


X 


10 2 


1.96 x 10 3 


1.92 


X 


10 3 


3.92 x 10 3 


65536 


7.82 


X 


10 2 


3.05 x 10 3 


2.96 


X 


10 3 


7.22 x 10 3 


131072 


1.45 


X 


10 3 




4.63 


X 


10 3 


1.34 x 10 4 


262144 


2.71 


X 


10 3 




7.36 


X 


10 :i 


2.50 x 10 4 


524288 


5.07 


X 


10 3 




1.19 


X 


10 4 


4.68 x 10 4 


1048576 


9.54 


X 


10 3 




1.97 


X 


10 4 


8.82 x 10 4 


2097152 


1.80 


X 


If) 4 




3.33 


X 


10 4 


1.67 x 10 5 



models is defined as the mass of stars with E > Edit, although 
the apocentre criterion is used in the simulations. The mass of stars 
satisfying the apocentre criterion is smaller than that of stars with 
E > .Ecrit, but shows a similar trend in time variation. 

Fig. [3] shows the half-mass time Tkaif, which is the time re- 
quired for a cluster to lose a half of its initial mass, as a function of 
the initial half-mass relaxation time t r h,i. Here the half-mass relax- 
ation time is defined by 



N^r 3 h /2 



t ^- - 138 - G i/2 m i/2 HjN) 
dSpitzerlll987t Chapt. 2) with 7 



(19) 

0.11 jGiersz&Heggie]|l994ah . 



The results are summarised also in Table Q] The FPf and TV-body 
models show good agreement over the whole range of TV where the 
comparison is made. The scaling Tiiaif x t?/ 4 gives a reasonable 
fit to the results of these models as Baumgardt (2001) found. 

The results of FPd models are also shown in Fig. [5] In 
these models the parameter v e = 2.5 is used for equation l[3} 
(Takahashi & Portegies Zwart 2000). The FPd models show clearly 
a different scaling from the other models; Thaif oc t r h expect for 
models with very short t^h (i- e - small TV ). 

In Fig. [4] FPf models with the energy criterion are compared 
with those with the apocentre criterion as well as the TV-body mod- 
els. We have set i/ c — 5 in the energy-criterion models so that 
their mass evolution reasonably agrees with that of the TV-body 
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the initial half-mass relaxation time t r h ; for the models shown in Fig. 




10 io 2 10 3 10 4 10 5 10 6 10 7 

Figure 5. Same as Fig.fJ] but FPf models with TV up to 2 30 are shown. The 
steeper dotted line represents the relation Th a if = t r h i. 

models for small TV. There is no significant difference between the 
energy-criterion models and the other models for t r h,i < 100, but 
the energy-criterion models tend to lose mass much faster as t r h,i 
increases. This indicates that the apocentre criterion is a better es- 
cape criterion for FPf models. 

As stated above, the results of the FPf models shown in Fig. [3] 
are reasonably well described by the scaling law Ti m if oc 
However we should not expect this scaling continues to hold in 
the limit of large TV. If this scaling continues, the half-mass time 
measured in the units of the half-mass relaxation time, Th a if/t r h,i, 
would go to zero as TV — > oo. This must be impossible because 
the mass loss is driven by two-body relaxation. In order to see the 
scaling of Thaif in the limit of large TV, we have calculated FPf 



Table 2. Half-mass times Ti m if given by FPf models for the initial condi- 
tions of Wo = 3 King models with very large TV. 





TV 


*rh,i 


Thaif 


Thaif Aril, i 


2 22 (« 


4.19 x 10 6 ) 


3.41 x 10 4 


5.77 x 10 4 


1.69 


2 23 (» 


8.39 x 10 6 ) 


6.47 x 10 4 


1.02 x 10 5 


1.57 


2 24 (as 


1.68 x 10 7 ) 


1.23 x 10 5 


1.82 x 10 5 


1.48 


2 25 (w 


3.36 x 10 7 ) 


2.35 x 10 5 


3.31 x 10 5 


1.41 


2 26 (w 


6.71 x 10 7 ) 


4.50 x 10 5 


6.09 x 10 5 


1.35 


2 27 (as 


1.34 x 10 8 ) 


8.62 x 10 5 


1.14 x 10 6 


1.32 


2 28 (ss 


2.68 x 10 8 ) 


1.65 x 10 6 


2.14 x 10 6 


1.30 


2 29 (« 


5.37 x 10 8 ) 


3.18 x 10 6 


4.07 x 10 6 


1.28 


2 30 (w 


1.07 x 10 9 ) 


6.12 x 10 6 


7.77 x 10 6 


1.27 



models with very large TV, TV = 2 22 r; 4.19 x 10 6 to 2 30 w 
1.07 x 10 9 , which are much larger than typical numbers of stars in 
globular clusters. The results of these models are shown in Table|2] 
and Fig. [5] In Fig. [5] we see that Thaif is nearly proportional to 
trh.i for very large TV clusters, say, for trh i > 10 5 or TV > 10 7 . This 
trend is more qualitatively shown in Fig.[6] where the change in the 
logarithmic slope, 

d log Thaif 

a =-Ti — , (20) 

dlog trh,i 

is plotted. The slope a approaches one as TV increases. The ratio 
Thaif/irh.i ~ 1.3 for our largest-TV models. 

We have performed simulations also for the initial conditions 
of Wo = 5 King models. The half-mass times of TV-body and FPf 
models for Wo = 5 are summarised in Table [3] and are plotted in 
Fig. [7] Here we find good agreement again. The same parameter 
v a = 7 is used for both the Wo = 3 and Wo = 5 clusters. In Fig. [7] 
the slope of the log t r h,i-log Thaif relation seems to be in between 
3/4 and 1. This point is further examined in subsection l3.3l 
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Figure 7. Same as Fig. [3] but for the initial conditions of Wo = 5 King 
models. FPf models with the apocentre criterion and ./V-body models are 
shown. 



Table 3. Half-mass times T[ la if given by TV-body and FPf models for the 
initial conditions of King models with Wo = 5. 



TV 


*rh,i 


Thaif 


Thaif 






(TV-body) 


(FPf) 


1024 


2.19 X 10 1 


3.89 x 10 2 


3.92 x 10 2 


2048 


3.82 x 10 1 


5.78 x 10 2 


6.07 x 10 2 


4096 


6.77 x 10 1 


9.51 x 10 2 


9.77 x 10 2 


8192 


1.22 x 10 2 


1.51 x 10 3 


1.61 x 10 3 


16384 


2.21 x 10 2 


2.54 x 10 3 


2.67 x 10 3 


32768 


4.04 x 10 2 


4.14 x 10 3 


4.49 x 10 3 


65536 


7.45 x 10 2 




7.62 x 10 3 


131072 


1.38 x 10 3 




1.31 x 10 4 


262144 


2.58 x 10 3 




2.28 x 10 4 


524288 


4.83 x 10 3 




4.03 x 10 4 


1048576 


9.09 x 10 3 




7.20 x 10 4 


2097152 


1.72 x 10 4 




1.30 x 10 5 




Figure 8. Same as Fig. [3] but FPf models with different functional forms 
of t e (E) oc (E - E clit )-P (J3 = 1, 2, 3) are compared. The dotted lines 
represent scalings t 2 / 3 , t^ 4 and t 4 ^ J , which are predicted by the simple 
steady-solution model for = 1, 2 and 3, respectively (see text). 

half-mass relaxation time for these runs as well as for the standard 
runs, where King models with Wo = 3 are used as initial condi- 
tions. The value of u c has been adjusted so that the non-standard 
models should have roughly the same half-mass times with those 
of the standard ones for lower TV; v Q — 7/3 and 7 x 3 for /3 = 1 
and 3, respectively. 

The results of the FPf models actually depend on /3, but the de- 
gree of the dependence is weaker than predicted by equation J2U . 
While this equation predicts the slopes 2/3, 3/4 and 4/5 for (9 = 1, 
2 and 3, respectively, linear least-squares fitting of the data in Fig.[8] 
gives the slopes 0.69, 0.72 and 0.75. When the fitting is done only 
for TV > 16384, the slopes are 0.75, 0.75 and 0.77. Thus the scal- 
ing law Ti la if oc t^ 4 is not a bad approximation in all the cases 
investigated here. This is not consistent with equation d2U . 



3.2 Dependence on the escape-time function 

iBaumgardtl ( 1200 ll) argued that the scaling Thaif oc can be ex- 
plained by a steady state solution of a simple model for the evo- 
lution of potential escapers (see equation (12) of his paper). His 
model adopts the escape time-scale t e of equation H3\ . If a differ- 
ent function is assumed for t e , his model predicts a different scaling 
law. It is shown that the scaling 

/3 + 1 

T ha i f OC t]g* (21) 

is obtained for t e oc (E — -Ecrit) - ^ (see Appendix A). It is inter- 
esting to see if this prediction is confirmed by the results of our FPf 
models. 

We have performed FP runs using a generalized form of equa- 
tion l fl5l , 

U{E) ~ ttid E C r\t ' 

with f3 = 1 and 3. Fig.[8]plots the half-mass time against the initial 



3.3 Dependence on the strength of the tidal field 

Tanikawa & Fukushig^ J2005h found that the dependence of Th a if 
on i r h,i is affected by the strength of the tidal field and that the 
logarithmic slope a, defined by equation l |20t . approaches unity as 
the strength of the tidal field decreases. In order to confirm their 
findings, we have calculated FPf models for the initial conditions 
where the initial tidal radius r t ,i is greater than the King cutoff 
radius tk (i.e. the radius at which the density drops to zero) for each 
value of Wo- On the other hand, all the models presented above are 
calculated for the initial conditions with r t ,i = nc. 

Table [4] lists the half-mass times for Wo = 3 King models 
with r t ,i/rK = 1.4, 2, 4 and 6, and Fig.[9]illustrates these results. 
In this figure the results for Wo = 3 and Wo = 5 King models 
with j"t,i/rK = 1 are also plotted. Note that the ratio tk(Wo = 
5)/ r K(Wo = 3) « 1.4. Fig.llOlshows the variation of a witht r h,i. 

The results shown in Figs. l9l and 1 101 confirm the findings of 
lTanikawa & Fukushige] fc005t) . The dependence of Ti m if on t^s 
does depend on the strength of the tidal field. In the limit of 
ft,i/n< — ► oo and TV — > oo, it is expected that a —5- 1. 
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Figure 9. Same as Fig. [5] but FPf models for the initial conditions of King 
models with r t ; > rx are compared with the cases of r t j = rjj. 




rh,i 



Figure 10. Logarithmic slope a = dlogTi la if/dlogt r h4 as a function of 
the initial half-mass relaxation time t r j, ; . The models are the same as those 
shown in Fig. [9] 



Note that the curve for Wo = 3 King models with rt.i/rK = 
1.4 lies very close to that for Wo = 5 King models with r t ,i/rK = 
1 in each of Figs. [9] and [To] This indicates that the mass-loss time- 
scale does not depend very much on the initial concentration of the 
cluster but is m ainly determined by the strength of the tidal field, 
as was found by Tanikawa & Fukushig 



3.4 Comparison with TV-body models: multi-mass clusters 

So far we have concentrated on single-mass clusters. Here we 
consider the evolution of multi-mass clusters comparing our FP 
models with the TV-body models of jGieles & Baumg ardtl J2008h . 



Table 5. Half-mass times T ha if given by TV-body {Gieles & Ba umgardl 
2008) and FPf models for the initial conditions of multi-mass King mod- 
els with Wo = 5 and r t ;/ric = 1. Three sets of the parameters (7, u c ) 
are used for the FPf models. 



TV 


Thalf 


^half 


Thalf 


Thalf 




(TV-body) 


(FPf) 


(FPf) 


(FPf) 






(0.11,7) 


(0.02, 7) 


(0.02, 40) 


1024 


1.14 x 10 2 


1.20 x 10 2 


1.69 x 10 2 


1.17 x 10 2 


2048 


1.74 x 10 2 


1.87 x 10 2 


2.54 x 10 2 


1.80 x 10 2 


4096 


2.69 x 10 2 


2.86 x 10 2 


3.75 x 10 2 


2.72 x 10 2 


8192 


4.35 x 10 2 


4.39 x 10 2 


5.59 x 10 2 


4.18 x 10 2 


16384 


6.70 x 10 2 


6.90 x 10 2 


8.57 x 10 2 


6.59 x 10 2 


32768 


1.06 x 10 3 


1.12 x 10 3 


1.36 x 10 3 


1.08 x 10 3 



10 c 



10* 



W =5, multi-mass 




FPf (7 = 0.1 =7) 
FPf (7 = 0.02,1/ =7) 
FPf (7 = 0.02,i/ e =40)_ 
N-body 

1 1 1 : 



10' 



10 4 

N 



10- 



Figure 11. Half-mass time T^aif as a function of the initial num- 
ber of stars TV for Wo = 5 King models with the IMF dN/dm oc 
m — 2.35 (m max /m m i n = 30). FPf models with three different sets 
of the parameters 7 and v e are compared with the TV-body models of 
iGieles & Baumgardtl f20o3) . 



They performed TV-body simulations of clusters on circular orbits 
around a point-mass galaxy. In their simulations the initial mass 
function (IMF) is given by dN/dm oc m -2 ' 35 with the ratio 
"imax/Wrain = 30. Stellar evolution is not considered in their 
simulations. The clusters initially have the density distribution of 
King models with Wo = 5. The ratio of the initial tidal radius to 
the King radi us rt.i/rK is varied from 1 to 8. The results of the 
simulations o f lGieles & Baumgardtl ( 2008) are summarised in their 
Table 1. Note that they use different notations from ours: rj is for 
the tidal (Jacobi) radius and rt is for the King radius. 

FP f models are calculated for the same initial conditions as 
those of IGieles & Baumgardtl d2008h . The results for rt.i/rx = 1 
are summarised in Tableland Fig.[TT] There the results of the FPf 
models with three differe nt sets of parameters 7 and v e are reported. 
lGiersz&Hegde1 ( ll994ah estimated the best value of 7 = 0.11 for 
single-mass clusters b y comparing TV-body m odels with FP and 
gas models. Similarly Giersz & Heggitj f 19961) obtained 7 = 0.02 
for multi-mass with the IMF dN/dm oc m _2 ' 5 (m max /'Timin = 
37.5). We have calculated FPf models for multi-mass clusters using 
these two values of 7. 
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Table 4. Half-mass times T| m if given by FPf models for the initial conditions of King models 
with Wq = 3 and r t j > rjc- 



N Thalf 7h a lf Thalf Thalf 





( r t,i/ r K 


= 1.4) 


( r t,i/ r K 


= 2) 


( r t,i/ r K 


= 4) 


( r t,i/ r K 


= 6) 


128 


1.62 x 


10 2 


2.71 x 


10 2 


6.86 x 


10 2 


1.11 x 


10 3 


256 


2.22 X 


10 2 


3.48 x 


10 2 


7.61 x 


10 2 


1.12 x 


10 3 


512 


3.02 x 


10 2 


4.52 x 


10 2 


8.95 x 


10 2 


1.25 x 


10 3 


1024 


4.44 X 


10 2 


6.48 x 


10 2 


1.24 x 


10 3 


1.73 x 


10 3 


2048 


6.93 x 


10 2 


1.01 x 


10 3 


1.97 x 


10 3 


2.78 x 


10 3 


4096 


1.12 x 


10 s 


1.68 x 


10 3 


3.33 x 


10 3 


4.83 x 


10 3 


8192 


1.87 x 


10 3 


2.84 x 


10 3 


5.82 x 


10 3 


8.64 x 


10 3 


16384 


3.14 x 


10 3 


4.89 x 


10 3 


1.02 x 


10 4 


1.54 x 


10 4 


32768 


5.33 x 


10 3 


8.49 x 


10 3 


1.80 x 


If) 4 


2.73 x 


If) 4 


65536 


9.16 x 


10 3 


1.49 x 


10 4 


3.16 x 


If) 4 


4.77 x 


If) 4 


131072 


1.59 x 


10 4 


2.64 x 


10 4 


5.56 x 


10 4 


8.34 x 


10 4 


262144 


2.80 x 


10 4 


4.72 x 


10 4 


9.89 x 


10 4 


1.48 x 


10 6 


524288 


4.99 x 


10 4 


8.54 x 


10 4 


1.78 x 


10 s 


2.64 x 


10 s 


1048576 


8.98 x 


10 4 


1.56 x 


10 5 


3.24 x 


10 s 


4.79 x 


10 s 


2097152 


1.63 x 


10 5 


2.87 x 


10 5 


6.02 x 


10 s 


8.99 x 


10 s 



Table 6. Half-mass times Tj, a if given by FPf models for the initial condi- 
tions of multi-mass King models with Wo = 5 and r^j/rK = 2, 4, 8. The 
adopted parameter set is (7, u c ) = (0.11, 7). 



N 


Thalf 


Thalf 


Thalf 




(r t ,i/r K = 2) 


(r t ,i/r K = 4) 


(r t ,i/r K = 8) 


1024 


3.34 x 10 2 


8.03 x 10 2 


1.81 x 10 3 


2048 


5.12 x 10 2 


1.17 x 10 3 


2.41 x 10 3 


4096 


7.79 x 10 2 


1.72 x 10 3 


3.35 x 10 3 


8192 


1.21 x 10 3 


2.65 x 10 3 


5.08 x 10 3 


16384 


1.95 x 10 3 


4.30 x 10 3 


8.33 x 10 3 


32768 


3.26 x 10 3 


7.32 x 10 3 


1.45 x 10 4 


Table 7. Same as Table[6] but the results of FPf models with thf 


set (7, 


u c ) = (0.02, 40) 


are listed. 




N 


Thalf 


Thalf 


Thalf 




(n,i/r K = 2) 


(r t ,i/r K = 4) 


(r t ,i/r K = 8) 


1024 


3.68 x 10 2 


9.28 x 10 2 


2.20 x 10 3 


2048 


5.53 x 10 2 


1.32 x 10 3 


2.87 x 10 3 


4096 


8.25 x 10 2 


1.89 x 10 3 


3.86 x 10 3 


8192 


1.26 x 10 3 


2.84 x 10 3 


5.63 x 10 3 


16384 


2.01 x 10 3 


4.53 x 10 3 


8.92 x 10 3 


32768 


3.55 x 10 3 


7.62 x 10 3 


1.52 x 10 4 




Figure 12 . Same as Fig. II II but FPf models are compared with the ./V-body 
models of Gieles & Baumgardt (2008) for the clusters with r t ; > r^. The 
adopted parameter sets for the FPf models are (7, v a ) = (0.11, 7) and 
(0.02,40). 



Fig. II II shows that the parameter set (7, f c ) = (0.11,7) 
adopted for single-mass clusters gives good fit to the iV-body mod- 
els also for multi-mass clusters. On the other hand the parameter 
set (7, u c ) — (0.02, 7) results in a clear deviation from the N- 
body models. If we stick to 7 = 0.02, the value of v c needs to be 
increased to about 40 in order to obtain good agreement with the 
TV-body models. We will discuss in more detail what values of the 
parameters we should choose in the next section. 

The results for the initial conditions wit h rt.i > 7~k are shown 
in Tab les [6] and [7] and Fig. [12] The results of iGieles & Baum gardt 
(2008) are not shown in these tables (see their Table 1). Fig. 1121 
shows that the FPf models with (7, v c ) = (0.11, 7) are in good 
agreement with the iV-body models for r t ,i/rK = 2 and 4. The 
FPf models with (7, v c ) = (0.02, 40) are a little farther to the 



iV-body models but still follow them rather well. However, for 
r t,i/Vi< = 8, a noticeable difference is observed between the FPf 
and iV-body models; in Fig. [12] the curve for the iV-body mod- 
els is approximately linear but the slopes of the curves for the FPf 
models apparently change with N . Neither parameter set repro- 
duces the results of the TV-body models as well as in the cases of 
rt,i/rK < 8. The reason for this discrepancy is not clear at present, 
but there is a possibility that very early core-collapse in the mod- 
els with r t ,i/r-K = 8 is, at least partially, responsible for it. The 
FPf model with (7, u c ) = (0.11, 7) and rt.i/Vic = 8 experiences 
core collapse (bounce) at t = 0.006Th a if for N — 1024, and at 
t = 0.03T ha if for N = 32768. The Coulomb logarithm may take 
different values for pre-collapse and post-collapse stages (see the 
next section), which affects the time-scale of the evolution of FP 
models. 
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4 DISCUSSION 

We have shown that FP models can well follow the mass evolu- 
tion of star clusters in a tidal field if a new scheme for treating 
potential escapers is implemented. This is the first time the effect 
of re-scattering of potential escapers has been taken into account 
in FP models. Although Takahashi & Portegies Zwart (1998, 2000) 
showed that anisotropic FP models are in good agreement with N- 
body models for the mass evolution of star clusters in a galaxy, the 
tidal field is treated as a tidal cutoff rather than an actual force field. 
In the present study we have found that our new FP models are in 
good agreement with iV-body models calculated with the inclusion 
of the tidal force field. Thus the new scheme has improved the ac- 
curacy of FP models. 

Baumgardt (2001) argued that some potential escapers are 
scattered back to lower energies before they leave the cluster and 
that this complicates the scaling of the mass-loss time. The suc- 
cess of our models is consistent with his argument. Actually our 
equation for potential escapers, equation J14b . can be regarded as a 
generalization of the equation of his toy model, his equation (12), 
used for explaining the scaling 7h a if oc ifh • 

The toy model o f lBaumgardtlfeOOll) is useful for giving us in- 
sight into the effect of potential escapers on the cluster evolution. 
On the other hand, the results presented in subsection 13. 2| have re- 
vealed the limitation of the model. When the energy dependence 
of the escape time is artificially changed from the true one, the toy 
model does not correctly explain the results of our FP models. This 
failure of the toy model is not a big surprise, because it is only a 
simplified model based on many assumptions, some of which are 
not very realistic. For example, our simulations show that an ex- 
act steady state is never established, but the toy model assumes a 
steady state. In addition, the scaling of the cluster lifetime depends 
on the strength of the tidal field, as found bv lTanikawa & Fukushigd 
d2005h and confirmed by the present study, but the toy model does 
not take account of the strength of the tidal field. 

Our FP models show good agreement with iV-body models 
not only for single-mass clusters but also for multi-mass clusters. 
However, we have encountered a difficulty in determining proper 
values of the two parameters, 7 and v c , in the FP models. As shown 
in subsection !3.4l the parameter set (7, v e ) = (0.11, 7) brings good 
agreement for both single-mass and multi-mass clusters. Since the 
escape time-scale t c given by equation dl5t is expected to be in- 
dependent of stellar mass, it is natural that the same value of the 
parameter u c is applicable to both single-mass and multi-mass clus- 
ters. 

On the other hand, th e valu e of 7 is expected to depend on 
the stellar mass function. iHenonl j 19751) argued theoretically that 
the value of 7 is generally smaller in multi-mass clusters than in 
single-mass clusters. B ased on the results of iV-body simulations, 
lGiersz&HeggieUl994al) obtained a value of 7 = 0.11 for isolated 
single-mass clusters, an dG lersz &Heggid l fl996l) obtained a much 
smaller value, 7 = 0.02, for isolated multi-mass clusters having an 
IMF similar to the IMF used in our simulations. 

When we adopt the value of 7 = 0.02 for multi-mass clusters, 
we have to use a much larger value of v c , p c = 40, than the best 
value of u e = 7 for single-mass clusters, in order to obtain good 
agreement with iV-body models. Thus we have not found a param- 
eter set satisfying both the independence of v c on the mass func- 
tion and the dependence of 7 on it. It needs further investigation 
to solve this incompatibility, but ev en the determination of 7 itself 
is not a simple task. For example, iGiersz & Heggiel Jl994bl) ob- 
tained the best value of 7 = 0.035 by examining the post-collapse 



evolution of iV-body models of isolated single-mass clusters. This 
value is much smaller than the value of 7 = 0.11 obtained for 
pre-collapse single-mass clusters. These results suggest that the 
value of 7 changes along with the e volution of clusters. It ma y also 
change with radius within a cluster dGiersz &Heggieill994j) . 

Fukushige & Heggie fcOOOl) theoretically estimated not only 
the energy dependence of the escape time-scale t e but also its nu- 
merical coefficient, which is given in their equation (9). If we ig- 
nore the difference between energy E and the lacobi integral Ej, 
their estimate for a Wo = 3 King model leads to a value of 
v c — 29. This is about four times larger t han our best value of 
u e = 7 for single-mass clusters. However, Fukushige & Heggie 
(2000) also did numerical experiments and found that their theoret- 
ical estimate of t c is too small; escape time-scales obtained from 
the numerical experiments are more than a few times larger than 
the theoretical on e. Therefore our value i/ e = 7 is not inconsistent 
with the result of iFukushige & H eggie ( 2000). On the other hand, 
our value of v c = 40 for multi-mass clusters with 7 = 0.02 is a 
little larger than their theoretical estimate. 

Another issue not addressed in the present paper is how the 
mass profile of the parent galaxy affects the results. In all the sim- 
ulations presented here we assume that th e parent galaxy is repre- 
sented by a point mass. On the other hand, Tanikawa & Fukushige 
d2010t) showed that the mass-loss time-scale depends on the mass 
profile of the parent galaxy; the time-scale increases as the mass 
profile gets shallower. Therefore we expect that the parameter i/ c 
depends on the mass profile of the parent galaxy. This issue will be 
examined in a future study. 



5 CONCLUSION 

In this paper we have developed new FP models of globular clus- 
ters in a steady galactic tidal field. Our FP models are novel in 
the method of treating escapers: potential escapers are allowed to 
experience gravitational scattering with other stars before they re- 
ally leave clusters. The new method has been devised in order to 
construct more realistic models of star clusters in a tidal field com- 
pared to simple tidal-cutoff models as in previous studies. The mass 
evolution of clusters in a tidal field does not simply scale with the 
relaxation time, and our FP models are in good agreement with N- 
body models in this respect. 

Our FP models include two parameters 7 and v e ; 7 is the 
numerical factor in the Coulomb logarithm ln(7iV) and v e ad- 
justs the speed of the tidal mass loss. We have determined the best 
values of v c for given values of 7 by comparing FP results with 
iV-body results. For single-mass clusters the best parameter set is 
(7, v e ) = (0.11, 7). This parameter set is applicable to multi-mass 
clusters as well, but another set (7, u c ) = (0.02, 40) does work 
equally well as long as multi-mass clusters are concerned. The pa- 
rameter v c is expected to depend on the mass profile of the parent 
galaxy, though a point-mass galaxy is assumed in all the simula- 
tions of the present paper. Further investigation is required for the 
determination of the best values of the parameters 7 and v c under 
various conditions. 

While FP models are generally thought to be less faithful mod- 
els of globular clusters than iV-body models, the present study has 
significantly improved the accuracy of FP models. An advantage of 
FP models is that they can be calculated much faster than iV-body 
models. Therefore FP models are particularly useful when we need 
to calculate a huge number of models. For example, when we try 
to specify the initial conditions of individual clusters, we have to 
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perform simulations for many sets of the initial conditions, because 
the parameter space to be searched is very large. We believe that 
our FP models is quite useful for such searching. 
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of escapers is nearly in equilibrium, equation dA2t shows that the 
width of the distribution is approximately given by 



AE 



and the number of escapers N c . 
is estimated to be 



(A3) 

NAE. The escape rate N eBC 



N n esc \ 
WrtJ 



fl+1 
■ 1 



t e (AE) tcsc 
Therefore the scaling of the half-mass time is given by 



Tha 



N 
N csc 



P + l 1 



(A4) 



(A5) 
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APPENDIX A: ESTIMATION OF THE SCALING OF THE 
CLUSTER LIFETIME 

We fo llow the arguments given by iBaumgardtl d200ll) and [Heggie 
d200lh in order to derive the scaling law of equation d — 1 1 > - 

Let E — (E — Edit) /l-Ecrit | and assume that the escape time- 
scale t c has energy-dependence such as 

t c (E) = t csc E~P (p > 0). (Al) 
Then Baumgardt's toy model is modified as 

9* = k L #n_ £ , p n_ 

Ot trf, dE 2 tcsc 

where n(E, t)dE is the number of stars with energies in the range 
(E, E+dE) and fci is a constant. If we assume that the distribution 



